5.3 构建环境因子背景矩阵

## 筛选环境变量的一些原则：
遵循以下原则：
参照Guisan[22]和Zhu[23]等的建议设定4条筛选原则来选择环境变量：
1、尽可能选择直接影响物种实际分布的近端作用变量和物种的生态需求密切相关的变量；2、基于物种生理耐受性的Jackniffe方法[24]，利用SDM(如Maxent)选择单因子建模更高解释度的变量[附录S2]。
3、为使建模结果更具有可解释性以及降低过拟合的可能性，需要去除变量间的强共线性关系(|r|>0.8)[25][附录S3]；
4、去除生物气候变量中部分温度和降雨结合的变量。

[23] Fan J Y , Zhao N X , Li M , et al. (2018)What are the best predictors for invasive potential of weeds? Transferability evaluations of model predictions based on diverse environmental data sets for Flaveria bidentis[J]. Weed Research.
[24] Steven J. P , Robert P.ER , Robert E. S. (2006)Maximum entropy modeling of species geographic distributions[J]. Ecological Modelling.

5.3.1 提取环境背景信息

library(raster)
library(dsimo)
## exs为数据矩阵，包含每个因子对应的点值信息；
> exs <- extract(predictors, bradypus)
summary(exs)

## 补充：
## 注意extract函数有多种用法，直接使用extract会导致采取的样点值为单元格的质心，可能会影响结果；
来自：
https://ecologicaconciencia.wordpress.com/
# 注意提取栅格值时，如果是在平原多水地区可使用线性平滑的方法进行提取，extract()有内部函数参数；
elev2<-resample(x=elev, y=ndvi, method="bilinear")

5.3.2 环境背景信息数据转换

## 可能会有时需要对分组数据做批量处理,如求回归,求极值,求均值等
## 使用aggregate()函数针对物种:list(sdmdata$species)
## aggregate(data,list=(),fun())
aggregate(iris[1:4],list(iris$Species),mean)

5.4 环境因子相关性分析

5.4.1 因子筛选主流方法

## 说明：
目前关于环境因子分析和筛选过程中，目前已经得到较认可的环境筛选方法是：
## 1、基于刀切法获取重要环境因子：
将所有环境因子纳入到建模分析钟，选择maxnet或者其他模型进行建模分析，根据建模的混淆矩阵评估建模结果的优劣；选择其中较好的建模结果(AUC.kappa)的结果模型提取主要贡献环境因子。
## 2、基于pearson或者spearman相关分析，绘制相关分析热图矩阵；
利用提取的环境因子信息矩阵，进行pearson分析，筛选环境变量相关性r>0.8以上的环境变量；
## 3、基于物种相关文献和背景环境，专家意见选择环境变量；
结合研究对象的生态分布信息，和植物生理生化需求选择合适的环境模糊因子；
## 4、比较和替换刀切法的环境因子，筛选建模因子
    因为刀切法中的建模环境因子仅基于数值统计，反映的是数据特征和数据波动的概率偏好，但不能反映的实际的环境因子需求；因此，使用3中筛选的模糊环境因子替换掉相关性较强的非相关生理因子进行建模；
    统计建模结果的auc或者其他统计变量，若模型仍有较低的遗漏率，则说明环境因子的替代是合理的；如果建模评价变差，则需要重新考虑和替换，原因可能有 如下几个方面：第一是专家意见判读不准确；第二是因子间关系非线性，有可能是高维相关；
第三是：模型在调参过程中选择的参数修正方式对线性关系不敏感；

## R 代码实现：
#  构建maxent基于全部环境变量的结果，查询建模结果的重要环境因子；
#  一般认为贡献率超过5%以上均可纳入模型中；
## 见 4.3构筑环境数据集合和5.1.5批量转asc文件；
## 见 6.1 sdm运行和sdm结果检验；
## 基于maxent到处结果，查询刀切法下的环境变量的重要性；
##3.2 查看刀切法构建的环境变量的敏感性(重要性)
plot(rattler.me)

## R 代码实现方法1：
# 构建pearson相关分析(pearson>0.75)：
## 去掉列名，仅保留各项环境因子：不需要进行数据标准化；
p1 <- cor(data1[,-1],method = "pearson")
library(pheatmap)  
## pheatmap：https://www.jianshu.com/p/1c55ea64ff3f
pheatmap(p1,cluster_cols = F,display_numbers = T,main = "环境因子相关性分析")
## 另外一种相关性组图的形式展示：
## 这种方法的缺点是当因子过多时,可视化效果变差;
pairs(sdmdata[,2:5], cex=0.1, fig=TRUE)
## 之后,基于上述两种分析的结果进行是实际的比较替换即可;


## R 代码实现方法2：
set.seed(123)
dd = as.data.frame(matrix(rnorm(1000),100,10))
# 计算相关系数及显著性：使用Hmisc包的rcorr()，可以得到correlation matrix的p值矩阵
library(Hmisc)#加载包D的
res2 <- rcorr(as.matrix(dd))
# 可视化
library(PerformanceAnalytics)#加载包
chart.Correlation(dd, histogram=TRUE, pch=19)

## 可视化的另外一种方法：弦图
## 哈哈哈，有空学习下；
https://www.sohu.com/a/400959520_120736615

5.4.2 因子筛选其他方法

## 利用方差膨胀因子筛选:
library(dismo)
predictors <- stack(files)
## 选择其中一个作为核心建模因子:
VIF.analysis <- ensemble.VIF(predictors, factors="xiz高程.tif ")
predictors

5.4.3 ENFA

## Biomapper是一套GIS和统计工具，旨在为任何种类的动物或植物建立栖息地适应性（HS）模型和地图。它以生态位生态因子分析（ENFA）为中心，该分析可以计算HS模型而无需缺勤数据。

## 软件为arcgis工具箱：
## https://www2.unil.ch/biomapper/index.html

刀切法、置换重要性及ENFA

## 因为需要构建四大洲的环境变量分析，因此这个工作需要在四大洲的层次上进行：
## load packages:
rm(list = ls())
setwd("C:\\Users\\admin\\Desktop")
##构建enm：
# install.packages("SDMtune")
library(SDMtune)
library(ENMeval)
library(raster)
library(rgdal)
library(maps)
library(mapdata)
library(dismo)
library(rJava)
library(maptools)
library(jsonlite)
library(glmnet)
library(maxnet)
library(rasterVis)
library(ggplot2)
##检查maxent是否可用：
dismo::maxent()

## 准备建模数据：
## 加载物种分布数据：####
## occs:

xh_sa <- read.csv("D:/XH/第二阶段/xh/na.csv")
head(xh_sa)
xh_sa <- xh_sa[,2:3]
####  加载环境数据：#####
## ENVS:
tifs <- list.files(path ="D:/XH/第三阶段_env/envsV3tif/",pattern="tif",full.names =T)
tiffs <- stack(tifs)


crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(tiffs) <- crs.geo  


## 根据物种分布数据构建建模范围 ########
## 以sa_xh为例：
occs <- xh_sa
xmin <- min(occs$longitude)
xmax <- max(occs$longitude)
ymin <- min(occs$latitude)
ymax <- max(occs$latitude)
bb <- matrix(c(xmin, xmin, xmax, xmax, xmin, ymin, ymax, ymax, ymin, ymin), ncol=2)
bgExt <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(bb)), 1)))
plot(env_bgExt)
# 增加buffer:0.5度~50km etc
bgExt <- rgeos::gBuffer(bgExt, width = 0.5)
## 利用环境数据掩膜：去除水面分布背景范围：
## 仅加载其中一张环境图层，用于裁剪-掩膜：
env_bgExt <- raster::mask(crop(tiffs$BIO17, bgExt),bgExt)
## 注意批量图层掩膜时，选择先用范围再裁剪的方法：
envs_bios <- crop(tiffs,env_bgExt)
## 生成用于建模的研究区域：
# plot(envs_bios$ARID)
# points(occs)
# plot(envs_bios$GDM10)
## 生成不包含研究区域的随机背景点：
bg.xy <- dismo::randomPoints(env_bgExt, 5000,p= occs)
bg_xy <- as.data.frame(bg.xy,colnames(c("long","lat")))

crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(envs_bios) <- crs.geo  
data <- prepareSWD(species = "xh_sa species", 
                   p =occs , a = bg_xy, 
                   env = envs_bios)
head(data@data)

## 去除无关变量：#########
## 方法1：maxnet：sdmtune：#####
## 评估PI JK
## 计算PI:
library(zeallot)  # For unpacking assignment
library(ENMeval)
c(train, test) %<-% trainValTest(data, test = 0.2, only_presence = TRUE, seed = 25)
maxnet_model <- train("Maxnet", data = train, fc = "lqp", reg = 1,iter=500)
auc(maxnet_model)


vi_maxnet <- varImp(maxnet_model, permut = 100)
vi_sa <- vi_maxnet 
## 计算：
## 
# the regularised gain of the variable when used
# in isolation, and the decrease in gain that occurs when
# a particular variable is omitted, was ranked
jk <- doJk(maxnet_model, metric = "auc", test = test)
jk_sa <- jk
## 方法2：BRT #########
library(dismo)
data(Anguilla_train)
head(Anguilla_train)
## 实际建模中需要提供的数据是：
## 对应研究区域分布点的点值提取结果；包含分布点和背景点，可以使用sduntune的框架来做；
data@pa
## 构建建模训练集：
## 提取发生数据：
sdmdata_sa <- cbind(data@data,data@pa)
dim(sdmdata_sa)
head(sdmdata_sa)
## 经检验，迭代到0.001时数据无法迭代成功；故选择0.01；
angaus.tc5.lr01 <- gbm.step(data=sdmdata_sa, gbm.x = 1:26, gbm.y = 27,
                            family = "bernoulli", tree.complexity = 5,
                            learning.rate = 0.01, bag.fraction = 0.5)

names(angaus.tc5.lr01)
## 查看变量的重要性：
brt_sa <- summary(angaus.tc5.lr01)
## 方法3：ENFA:  ##########
library(CENFA)
occ_sa <- data@coords[data@pa==1,]
head(occ_sa)
## 点转shp，空间数据库形式：
head(xh_sa)
library(sp)
coordinates(xh_sa) <- ~ longitude +latitude ##这里的lon和lat是清理数据集中指定的经纬对应列
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(xh_sa) <- crs.geo  

## 需要构架分布点的最小包络：
### 构建mcp：
## install.packages("adehabitatHR")
library(adehabitatHR)
##结果证明，使用kde，确实比mcp好很多呀；
occ_kde <- kernelUD(xh_sa, h="LSCV", hlim = c(0.5, 1.5))
occ_kde_poly <- getverticeshr(occ_kde, percent = 90) 

plot(turtle.kernel.poly, col = occ_kde$ud)
points(xh_sa)

mod.enfa <- enfa(x = envs_bios, s.dat = occ_kde_poly)
enfa_sa <- mod.enfa
glc <- GLcenfa(x = envs_bios)
scatter(x = mod.enfa, y = glc)
mod.cnfa <- cnfa(x = envs_bios, s.dat = occ_kde_poly)
cnfa_sa <- mod.cnfa

sink("./sa.txt")
vi_sa
jk_sa
brt_sa
cnfa_sa
sink()

其他筛选环境变量的方法

## 筛选变量的一种方法：
https://academic.oup.com/jmammal/article/96/3/511/905515
# 在这篇文章中提供一个思路给定协方差，通过评估不同协方差组合所构建的AIC指数，以及不同AIC指数与步长之间的交互关系来最终评估协方差在参与建模中的重要性；